A two-species continuum model for aeolian sand ripples 
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O , We formulate a continuum model for aeolian sand ripples consisting of two species of grains: a 

^^ ' lower layer of relatively immobile clusters, with an upper layer of highly mobile grains moving on 

top. We predict analytically the ripple wavelength, initial ripple growth rate and threshold saltation 
flux for ripple formation. Numerical simulations show the evolution of realistic ripple profiles from 
initial surface roughness via ripple growth and merger. 
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Aeolian sand ripples are formed by the action of the 
wind on the sand bed in the desert and at the seashore. 
They have also recently been observed on Mars ||l| . Aeo- 
^ . lian ripples are a few centimetres in wavelength and their 
"t^ I crests lie perpendicular to the prevailing wind direction. 
r^ . Bagnold (1941) made an influential early study M. He 
I ' identified the importance of saltation, where sand grains 
■ O ■ are entrained by the wind, and whipped along the sand 
C ] bed, colliding with it at high speed (of order < Im/s [^), 
and causing other grains to jump out of and along the 
bed, thus sculpting ripples. The impact angle remains 
roughly constant at about 10° — 16° despite the gusting 
K^ ' of the wind. The stoss slope of the ripple lies in the range 
fvi . 8°-10° [Q. The lee slope is composed of a short straight 
l/~) ' section near the crest at an angle of about 30° — 34° to 
C^l . the horizontal fj] 0, followed by a longer and shal- 
^^ ■ lower concave section. The deposition of grains on the lee 
slope of sand dunes leads to oversteepening and avalanch- 
iiig lo 101' which maintains the lee slope at an angle of 
around 30° — 34° near the crest. A similar, if less dra- 
d ' matic, mechanism is likely to hold for ripples. 
S . Numerical simulations of ripples and dunes based on 

I , ' tracking individual sand grains or on cellular automata 
^5 , have in the past few years yielded good qualitative agree- 
Q ' ment with observations [g[. However, these methods are 
Q . computationally expensive. Continuum models provide 
l^ ' a complementary approach allowing faster calculation of 
• i-H , ripple evolution, and the possibility of obtaining certain 
^\ information, such as the preferred ripple wavelength, an- 
alytically. Anderson [pi produced an analytical model of 
' the initial generation of ripples from a flat bed. A one- 
species analytical continuum model was formulated by 
Hoyle and Woods ||l^]. In the current Letter we extend 
the one-species model [|l3] to include relaxation effects, 
inspired by models of sandpiles [ pT| , in order to obtain 
more realistic ripple profiles and to predict the ripple 
wavelength and speed. As our work was nearing com- 
pletion we learned of other work in progress on ripple 
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formation [|2| similar in spirit to our own, but with im- 
portant differences in implementation. 

We consider two-dimensional sand ripples: following 
earlier work on sandpiles Q Q we assume that the 
effective surface of the ripple comprises a 'bare' surface 
defined by the local heights of clusters h{x,t) sheathed by 
a thin layer of flowing mobile grains whose local density is 
p{x, t), with X the horizontal coordinate, and t time. The 
ripple evolves under the influence of two distinct types 
of mechanism. Firstly, the impact of a constant flux of 
saltating grains, knocks grains out of the 'bare' surface, 
causing them to hop along the ripple surface and land in 
the layer of flowing grains. This is the underlying cause 
of ripple formation. Secondly, the ripples are subject to 
intracluster and intercluster granular relaxation mecha- 
nisms which result in a smoothing of the surface. 

Following jl^l we consider grains to be bounced out of 
the 'bare' surface by a constant incoming saltation flux, 
which impacts the sand bed at an angle (3 to the horizon- 
tal. These hopping or 'reptating' [|l^ grains subsequently 
land in the flowing layer. The saltating grains are highly 
energetic and continue in saltation upon rebounding from 
the sand bed. We assume that the number N{x, t) of sand 
grains ejected per unit time, per unit length of the sand 
bed, is proportional to the component of the saltation 
flux perpendicular to the sand surface, giving 

,x. X T . / ^N Jfsin /3 -f /la; COS /3) 
N{x,t) ^ J sm{a + (3)= ^ (^_^ ^2)1/2 ' (1) 

where a = tan'^{hx)j and J is a positive constant of 
proportionality. We assume that each sand grain ejected 
from the surface hops a horizontal distance a, with prob- 
ability p{a), and then lands in the flowing layer. We con- 
sider the flight of each sand grain to take place instanta- 
neously, since ripples evolve on a much slower timescale 
than that of a hop. The hop length distribution p{a) can 
be measured experimentally Q [^ , so we consider it to 
be given empirically. It is possible that p{a) and hence 
the mean and variance of hop lengths, could depend upon 
factors such as wind speed. The number Sno{x, t) of sand 
grains leaving the surface between positions x and x + 5x 
in time (5i, where 5x and 5t are infinitesimal is given 
by 5no{x^ t) — N{x, t)Sx6t. The change 6h in the surface 
height satisfies Sx6h{x, t) = —ap5no{x, i), where Op is the 
average cross-sectional area of a sand grain. In the limit 
(5t — > we find that the contribution to the evolution 
equation for h{x,t) from hopping alone is 

, ,w ,. sin/3-Hft,a;C0s/? 

ht = ^apN{x, t) = -QpJ (i + fe2)i/2 ■ (2) 

There may be regions on the sand bed which are shielded 
from the incoming saltation flux by higher relief upwind. 
In these regions there will be no grains bounced out of 
the surface, and there will be no contribution to the ht 
equation from hopping. The number 5ni{x,t) of sand 



grains arriving on the layer of flowing grains between 
positions x and x + 5x in time 5t is given by 

/+CX; 
p{a)N{x — a, t)daSx5t. (3) 

-C30 

The change in depth of the flowing layer satisfies 
SxSp{x,t) = apSni{x,t), and hence the contribution to 
the evolution equation for the flowing layer depth from 
hopping alone is 

/■+°° sin/3 + /i:r(a;-a,t)cos/3 

We incorporate diffusive motion jl^l as well as pro- 
cesses governing the transfer between flowing grains and 
clusters, following ^^, leading to the equations 

sin P + hx cos /3 
ht= Dhh^^~T(x,t)~ apJ——ji^yj:^, (5) 

pt= Dpp^x + x{pK)x + T{x, t) 

J f^°° I ^smP + hj;{x-a,t)cosf] 

where Z3/,,, Z?p and x ^^^ positive constants and where 
r(x, i), which represents the transfer terms, is given by 

T{x,t) ^ -Kphxx + >^pi\hx\ -tana), (7) 

for < I /la; I < tana and by 

rr/ ^\ u , I'd /la; I -tana) 

"^("'^) = '"^"-+ (tan^^/ig)V2 ' (8) 

for tana < |/ia;| < tan7, with k, A and i^ also positive 
constants. The term Dhh^x represents the diffusive re- 
arrangement of clusters while the term Dpp^x represents 
the diffusion of the flowing grains. The flux-divergence 
term x{phx)x models the flow of surface grains under 
gravity. The current of grains is assumed proportional to 
the number of flowing grains and to their velocity, which 
in turn is proportional to the local slope to leading or- 
der |10|] . The —Kphxx term represents the inertial filling 
in of dips and knocking out of bumps on the 'bare' sur- 
face caused by rolling grains flowing over the top. The 
\p{\hx\ — tana) term represents the tendency of flowing 
grains to stick onto the ripple surface at small slopes; it is 
meant to model the accumulation of slowly flowing grains 
at an obstacle. Clearly for this to happen the obstacle 
must be stable, or else it would be knocked off by the on- 
coming grains, so that this term only comes into play for 
slopes less than tan a, where a is the angle of repose. The 
term z/(|/ia;| — tana)(tan^ 7 — h'^)~^^'^ represents tilt and 
avalanching; it comes into play only for slopes greater 
than tan a and models the tendency of erstwhile stable 
clusters to shed grains into the flowing layer when tilted. 
This shedding of grains starts when the surface slope ex- 
ceeds the angle of repose. For slopes approaching the an- 
gle 7, which is the maximum angle of stability, the rate 



of tilting out of grains becomes very large: an avalanche 
occurs. This term, among other things, is a novel repre- 
sentation of the well-known phenomena of bistability and 
avalanching at the angle of repose IItI . 

We renormalise the model equations setting x — > 
xox, t — > tot, a —> xqu, p -^ Pop, h — > hoh, 
where xo = Dh/apJ cosP, to = Dh/{apJ cos PY , ho — 
Dhtan j/upJ cos (3, po — UpJ sinP/Xtana, which gives 
for < hx < tan a/ tan 7 

/, . N, tan/3 /,, , tanaX ^, , ,„, 
ht^ (1 + kp)h,., - p—^ [\hx\- - fix), (9) 



tan a \ tan 7 

ho ( ^ , tan/? / tana 

Pt= — < -uphxx + p- \hx 



Po [^ tana \ tan 7 



p{a)f{x - a)da + ^Pxx + x{phx)x, (10) 

Po J-00 ^h 

and for tan a/ tan j < hx < 1 

, n , - \u ft \ ^(l^:r| -tan a/ tan 7) 

/iQ f ^ i> (I /la; I — tan a/ tan 7) 

/3t= — { -uphxx + 



Pol ' {^-hiyi'^ 

H — ^/ p{a)f{x-a)da+ —f-pxx + xipfT'x)x, (12) 

Po J-00 ^?i 

where the tildes have been dropped and where 

f{x)^(hx + '£^yi + hltan'j)-'/^ (13) 

and K = npo/Dh, v = vto/ho, x = X^o/Dh- Wher- 
ever the sand bed is shielded from the saltation flux, the 
hopping term must be suppressed by removing the term 
—f{x) in the ht equation. 

Close to onset of the instability that gives rise to sand 
ripples, the slopes of the sand bed will be small, since 
surface roughness is of small amplitude; hence the regime 
< tan a/ tan 7 is relevant. There are no shielded regions 
at early times, since the slope of the bed does not exceed 
tan/3. Note that hx = 0, p = 1 is a, stationary solution 
of equations ( §) and ( |l^). Setting h = he'^^+i^^ and 
p = 1 -f pe'^t+'^kx ^ where /i ^ 1 and p -^ 1 are constants, 
linearising, and Taylor-expanding the integrand gives a 
dispersion relation for a in terms of k. The presence 
of the I hx I term means that strictly we are considering 
different solutions for sections of the ripple where hx > 
and sections where hx < 0, but in fact the effect of the 
\hx\ terms appears only as a contribution ± tan /3/ tan a 
to the bracket (1 ± tan/3/tana), where the hx > case 
takes the -t- sign and hx < case takes the — sign. Since 
typically we have tan /3/ tan a ^ 1, the two solutions 
will not be very different. One growth rate eigenvalue is 
given by cr = — /iq tan/3/po tan7 + 0{k) and is the rate 
of relaxation of p to its equilibrium value of 1. To 0{k'^) 
the other eigenvalue is 



cr = (a - 1 - xPQ/ho)k'^ + iAk^ + Bk'^ (14) 

where 

A^-^^+^'-^Jl + X^--a-^)(l±'-^], (15) 

2 /iQ tan p \ ho Dh J \ tan a J 

S = _l^_l^^!^fl±^„^^f,_l„;e? + #^')|a + ^+?^fl±^'|, (16) 
6 2 fto tan p \ tan a J iiq tan p \ h^ Dh J \ /lo tan p \ tan a J I 

and where (.) denotes /_ {.)p{a)da. We have neglected 
higher order terms as we are looking for long wave modes 
where |fc| is small, since short waves are damped by the 
diffusion terms. Sand ripples grow if a > 1 + xPo/^o, 
which is equivalent to requiring that aap J cos (3 > D^ + 
Xflp Jsin/3/A tana holds in physical variables, giving a 
threshold saltation flux intensity for ripple growth. This 
is in agreement with the threshold found in [Q. Since 
B is negative (/3 < a), the fastest growing mode has 
wavenumber k^ = —(a — 1 — xpQ/ho)/2B with growth 
rate a = —{a — 1 — xPo/ho)^/4:B. The allowed band of 
wavenumbers for growing modes is < A;^ < — (a — 1 — 
XPo/ho)/B. The wave speed is given by c = —Ak"^ > 0; 
it is higher for larger fc^ which implies that shorter waves 
move faster, as indeed was seen in the numerical simu- 
lations described below. The speed is higher for h^ > 
than for hx < 0, leading to wave steepening. 



The renormalised model equations ( W)-{ 12) were in- 
tegrated numerically using compact finite differences pq ] 
with periodic boundary conditions. The —f{x) term 
in the ht equation was suppressed in shielded regions. 
We used a normal distribution for the hop lengths with 
mean a and variance s^. In the run illustrated, we chose 
a = 3.1, s = 0.1, Dp/Dh - 1.0, /iq/po = 20.0, x = 0.1, 
i> = 1.0, k = 0.1, /3 = 10°, a = 30° and 7 = 35°. 
The angles were chosen to agree with observational ev- 
idence, the ratio /lo/po to ensure a thin layer of flow- 
ing grains, and the remaining parameters to allow ripple 
growth. The output was rescaled back into physical vari- 
ables using Dh — 1.0 and A — 10.0. The initial condi- 
tions for the dimensionless variables were h — l.Q + O.lrjh, 
p = 0.95 -|- 0.1?7p, where r^h and rjp represent random 
noise generated by random variables on [0, 1) in order 
to model surface roughness. In this case B = —8.47 and 
A = —5.26 (taking the minus sign in the brackets), giving 
a preferred wavenumber of fc = 0.352, and a wave speed 
of c = 5.26fc^. The length of the integration domain was 
chosen to be ten times the linearly preferred wavelength. 

Figure H shows the surface height at time t = lO.OAi, 
where Ai = 2.78. Note the emergence of a preferred 
wavelength, with wavenumber k k. 0.457 lying in the 
permitted band for growing modes predicted by the lin- 
ear stability analysis and arrived at by a process of ripple 
merger. The wave speed close to onset was also measured 
and found to be c « (4.13 ± 0.39)fc^, which is reasonably 
close to the predicted value. Ripple merger typically oc- 
curs when a small fast ripple catches up and merges with 



a larger slower ripple (figure ||) , the leading ripple trans- 
ferring sand to its pursuer until only the pursuer remains. 
Occasionally a small ripple emerges from the front of the 
new merged ripple and runs off ahead. Figure shows the 
surface height h at time t — 89Ai, with one shallow and 
one fully developed ripple. Note the long shallow stoss 
slopes, and the shorter steeper lee slopes with straight 
sections near the crests and concave tails. The leftmost 
ripple has a maximum stoss slope angle of 3.3°, and a 
maximum lee slope of 9.2°, whereas the more fully devel- 
oped rightmost ripple has a maximum stoss slope angle 
of 24.8° and a maximum lee slope angle of 33.5°, which 
lies between the angle of repose and the maximum angle 
of stability. The height to length ratio of the ripples is 
in the range 1:8 - 1:22, which is in reasonable agreement 
with observations [y. In the long time limit, we would 
expect sand ripples to grow until the maximum lee slope 
angle reaches an angle close to tan 7. In reality, there 
is only a relatively shallow layer of loose sand available 
for incorportion into ripples, and this together with the 
maximum slope condition will determine the size of the 
fully-developed ripples. 

In summary, we have formulated an analytical con- 
tinuum model for aeolian sand ripples using a two- 
species model embodying intracluster and intercluster re- 
laxation, in a description that leads naturally to bistable 
behaviour at the angle of repose, and its cutoff at the 
angle of maximal stability m^ . We have predicted ana- 
lytically the preferred ripple wavelength, the wave speed 
and the threshold saltation flux required for ripples to 
form. Our numerical simulations show the development 
of realistic ripple profiles from initial surface roughness 
via growth and ripple merger. 
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FIG. 1. The surface height h at time t = lO.OAt. 
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FIG. 2. A sequence of profiles showing a small ripple catching up and merging with a larger npple. Sand is transferred from 
the larger to the smaller ripple until only the latter remains. The profiles are shown at times t — ll.OAi (solid line), t = 12.0At 
(dashed line) and t = 13.0At (dot-dash line). The later profiles are each offset by one additional unit in height. 
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FIG. 3. The surface height h at time t 
close to the crests. 



89At, showing fully developed ripples. Note the straight segments on the lee slopes 



